--- title: Basics keywords: fastai sidebar: home_sidebar nb_path: "nbs/01_basics.ipynb" ---
{% raw %}
{% endraw %}

Fastai 2 provides support to read and display medical images using pydicom and pillow, however only 2D images can be read. Faimed3d provides tool to load an preocess 3d medical data

Create a DICOM base class

Fastai already provides PILDicom and TensorDicom. Since these objects are designed to handle 2D data, using then might lead to problems in later tasks. So a custom Tensor class will be created to be able to use some of the very handy fastai functions.

freqhist_bins and hist_scaled are needed to for image normalization

The retain_type function from fast ai can be used, to retain the TensorDicom3D type.

retain_type(t, typ = TensorDicom3D)

Read medical 3D images

3D data is stored in a variety of data formats (DICOM, NIfTI, NRRD, Analyze, ...) which should be supported by the image loader. Also DICOM data is often stored as individual slices (DICOM series) and not as a volume. SimpleITK is a powerfull library which can handle many data formats, including all of the above mentioned and many more. It will therefore be used to read the 3D volumes.

{% raw %}

readdcm_3d[source]

readdcm_3d(fn:Path'>, <class 'str'>), div=None, return_scaled=False, return_normalized=False)

Opens a 3D medical file or DICOM series and returns a tensor Args: fn: file name as string or pathlib.Path div: new max for rescaling pixel values. If None, pixels will be rescaled between 0 and 1.

{% endraw %} {% raw %}
{% endraw %}

Display 3D images

Data visualisation is very important. Analogous to a radiology workstation it should be possible to display the images in axial, coronal or sagittal reformation.

{% raw %}

show_image_3d[source]

show_image_3d(t:Tensor'>), axis:int=0, figsize:int=(15, 15), cmap:str='bone', nrow:int=10, alpha=1.0, return_grid=False, add_to_existing=False, **kwargs)

Plots 2D slices of a 3D image alongside a prior specified axis. Args: t: a 3D numpy.ndarray or torch.Tensor axis: axis to split 3D array to 2D images figsize, cmap: passed to plt.imshow nrow: passed to torchvision.utils.make_grid return_grid: Whether the grid should be returned for further processing or if the plot should be displayed. add_to_existing: if set to true, no new figure will be created. Used for mask overlays

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
im = readdcm_3d('/media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM')
im.shape, show_image_3d(im, axis = 0, nrow = 10)

((31, 736, 736), None)
{% endraw %}

Sometimes multiple 3D images (e.g. a batch) need to be displayed. With a wrapper for show_image_3d this is conveniently possible.

{% raw %}

show_images_3d[source]

show_images_3d(t:Tensor, axis:int=0, figsize:int=(15, 15), cmap:str='bone', nrow:int=10, alpha=1.0, return_grid=False, add_to_existing=False, **kwargs)

displays multiple 3D images (e.g. a batch) by flattening the 4th dimension of the tensor and then calling show_image_3d

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
im.shape

(31, 736, 736)
{% endraw %} {% raw %} 
show_images_3d(torch.stack((im,im)), axis = 0, nrow = 31, figsize = (25, 15))
{% endraw %}

Crop and resize 3D images

Often the important areas of the images are also located centrally. For example in CT scans, the skin should always be included in the image, which means that the outer pixels will display only air or the scanner table. These areas can be cropped, effectivly reducing the size of the image without reducing the resolution.

{% raw %}

crop_3d[source]

crop_3d(t:Tensor, margins:(<class 'int'>, <class 'float'>), perc_margins=False)

Often the regions of interest in the DICOM/NIfTI are centrally locaed in the image. Cropping the margins substantially reduces image size, which makes it possible to keep a higher resolution in the relevant regions. However, absolute cropping of pixel margins might be a problem with differen image resolutions for the same sequenze (e.g. T2 images can be 736x736 px or 448x448). ropping the same margin for differen pixel resolutions could result in unwandet zooming effects. Cropping a percentage margin might be more beneficial.

Args: t = torch Tensor object margins = a tuple of pixels to crop at both sides of each dimension. Either only one value for each dimension is given for a symmetrical crop or two values are given for asymmetric crop. perc_margins = whether the margins should be treated as absolute values or precentage values

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
cropped_im_abs = crop_3d(im, ((0,0),(100,100),(100,100)), perc_margins=False)
cropped_im_perc = crop_3d(im, ((0.,0.),(0.1,0.1),(0.1,0.1)), perc_margins=True)

im.shape, cropped_im_abs.shape, cropped_im_perc.shape

((31, 736, 736), (31, 536, 536), (31, 590, 590))
{% endraw %} {% raw %} 
show_image_3d(cropped_im_abs)
{% endraw %}

Medical data has different resolutions. Most CT scans will be at 512x512 px but for MRI the resolution varies. Also the number of slices is different for every scan. To train a neural network all input needs to be the same size, thus a function to resize the images alongside each axis is needed.

{% raw %} 
# works well but is overly complicated. Using F.interpolate is much simpler and is implemented as a tfms. 
def resize_3d_tensor(t: Tensor, new_shape: int):
    """
    A function to resize a 3D image using torch.nn.functional.grid_sample
        
    Args:
        t (Tensor): a Rank 3 Tensor to be resized
        new_dim (int): a tuple with the new x,y,z dimensions of the tensor after resize
        
    Returns:
        a resized `Tensor`. Note that the header with metadata is lost in this process. 
        
    Procedure: 
        [fake RGB] -> [batch dim] -> [flow field] -> [resampling] -> [squeeze]
        
        fake RGB: 
            Create a fake RGB 3D image through generating fake color channels.
        btach dim: 
            add a 5th batch dimension
        flow field: 
            create a flow-field for rescaling:
                a. create a 1D tensor giving a linear progression from -1 to 1
                b. creat a mesh-grid (the flow field) from x,y,z tensors from (a)
                
                Taken form the offical Pytroch documention: 
                    Given an input and a flow-field grid, computes the output using input values and pixel locations from grid.
                    In the spatial (4-D) case, for input with shape (N,C,Hin,Win) and with grid in shape (N, Hout, Wout, 2), the output will have shape (N, C, Hout,Wout)

                    In the case of 5D inputs, grid[n, d, h, w] specifies the x, y, z pixel locations for interpolating output[n, :, d, h, w]. 
                    mode argument specifies nearest or bilinear interpolation method to sample the input pixels.
                
        resampling:  
            resample the input tensor according to the flow field
        squeeze: 
            remove fake color channels and batch dim, returning only the 3D tensor  
    """
    
    if type(new_shape) is tuple and len(new_shape) == 3:
        z,x,y = new_shape # for a reason, I do currently not understand, order of the axis changes from resampling. flipping the order of x,y,z is the current workaround
    else:
        raise ValueError('"new_dim" must be a tuple with length 3, specifying the new (x,y,z) dimensions of the 3D tensor')
    
    t = t.unsqueeze(0) # create fake color channel
    t = t.unsqueeze(0).float() # create batch dim    

    x = torch.linspace(-1, 1, x) # create resampling 'directions' for pixels in each axis
    y = torch.linspace(-1, 1, y)
    z = torch.linspace(-1, 1, z)

    meshx, meshy, meshz = torch.meshgrid((x, y, z)) # 
    grid = torch.stack((meshy, meshx , meshz), 3) # create flow field. x and y need to be switched as otherwise the images are rotated. 
    grid = grid.unsqueeze(0) # add batch dim
    out = F.grid_sample(t, grid, align_corners=True, mode = 'bilinear') # rescale the 5D tensor
    return out.squeeze().permute(2,0,1).contiguous() # remove fake color channels and batch dim, reorder the image (the Z axis has moved to the back...)
{% endraw %} {% raw %} 
resized_im = resize_3d_tensor(im, (100,50,50))
resized_im.shape, show_image_3d(resized_im)

((100, 50, 50), None)
{% endraw %}

Normalize images

Especially for MRI the pixel values can vary between scanner types. However, the imagenet stats, as provided by fastai should probably not be used, as they are specfic for color images and not MRI images. The optimal solution would probably be to calculate the stats on the present dataset.

Statistics of intensity values, pooled across the whole training dataset (raw, unscaled pixel values and pixel values scaled between 0 and 1)

Statistic ADC T2 T1 map
Pooled Maximum of intensity values 3064.2500 1369.5652 4095
Pooled Minimum of intensity values 0. 0. 0.
Pooled Mean of intensity values 511.1060 259.1454 740.8268
Pooled SD of intensity values 488.8707 190.2448 688.8238
Scaled Maximum of intensity values 1. 1. 1.
Scaled Minimum of intensity values 0. 0. 0.
Scaled Mean of intensity values 0.1675 0.1918 0.1809
Scaled SD of intensity values 0.1599 0.1409 0.1682

However, just scaling wiht one mean and std might not be the optimal solution (see the excelent Kaggle kernel form Jeremy Howard why). Although for MRI images it might be ok, as the pixel values are closer together.

{% raw %} 
im1 = resize_3d_tensor(readdcm_3d('/media/ScaleOut/prostata/data/dcm/A0042280702/T2/DICOM'), (10, 100, 100))
im2 = resize_3d_tensor(readdcm_3d('/media/ScaleOut/prostata/data/dcm/A0041912240/T2/DICOM'), (10, 100, 100))
im1.shape, im2.shape

((10, 100, 100), (10, 100, 100))
{% endraw %} {% raw %}

class TensorDicom3D[source]

TensorDicom3D(*args, **kwargs) :: Tensor

Base class for 3D Dicom Tensor. Inherits from torch.Tensor

{% endraw %} {% raw %}

load_image_3d[source]

load_image_3d(fn:Path'>, <class 'str'>))

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
im = TensorDicom3D.create('/media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM')
{% endraw %} {% raw %} 
im.metadata['table']
tags /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0030.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0029.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0028.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0027.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0026.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0025.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0024.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0023.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0022.dcm ... /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0009.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0008.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0007.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0006.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0005.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0004.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0003.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0002.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0001.dcm /media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM/0000.dcm
0 0008|0005 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ... ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100 ISO_IR 100
1 0008|0008 ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ... ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT ORIGINAL\PRIMARY\M\NORM\DIS2D\SH\FIL\MFSPLIT
2 0008|0012 20181205 20181205 20181205 20181205 20181205 20181205 20181205 20181205 20181205 ... 20181205 20181205 20181205 20181205 20181205 20181205 20181205 20181205 20181205 20181205
3 0008|0013 201257.045534 201256.875414 201256.721304 201256.560189 201256.394070 201256.231956 201256.069840 201255.909727 201255.739605 ... 201254.137467 201254.029389 201253.915308 201253.788218 201253.668132 201253.541042 201253.412951 201253.284860 201253.132752 201253.012667
4 0008|0016 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 ... 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4 1.2.840.10008.5.1.4.1.1.4
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
102 0040|0245 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000 ... 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000 155756.000000
103 0040|0253 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 ... SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942 SIb5f61ffe5db942
104 0040|0254 MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken ... MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken MRT Becken
105 0040|2017 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000 ... 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000 050A008651890000
106 2050|0020 IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY ... IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY IDENTITY

107 rows × 32 columns

{% endraw %}

Segmentation Masks need an own class, because image transforms can only be applied to prior specified classes. Thus, transforms altering the pixel value (e.g. add noise, change brightness) will not be applied to Masks if the class differs from the original image.

{% raw %} 
class TensorMask3D(TensorDicom3D):
    "Base class for 3d Segmentation, inherits from TensorDicom3D"
    @classmethod
    def create(cls, fn:(Path,str,Tensor,ndarray), **kwargs)->None:
        "open the mask, keep as float"
        if isinstance(fn,ndarray): return cls(fn)
        if isinstance(fn, Tensor): return cls(fn)
                    
        # nifti files differ in pixel orientation from DICOM.
        # This might be due to differnt pixel orientations. 
        # currently this problem was solved otherwise, so the create function does not differ from TensorDicom3D. 
        # However future flipping and Rotation might be implemented here. 
            
        return cls(load_image_3d(fn))
{% endraw %} {% raw %} 
px1 = TensorDicom3D(im1.flatten())
plt.hist(px1.numpy(), bins = 40)
plt.hist(px1.hist_scaled().numpy(), bins = 40)
plt.title("Distribution of pixel values in T2 images original and scaled") 
plt.show()
{% endraw %}

The pixels are now all nearly equally distributed.

{% raw %} 
im1 = TensorDicom3D(im1)
im2 = TensorDicom3D(im2)
{% endraw %} {% raw %} 
show_images_3d(torch.stack((im1, im2, im1.hist_scaled(), im2.hist_scaled()), axis = 0), nrow = 10)
{% endraw %}

I belive both steps should executed after reading the images, even before cropping and resizing, thus they will be directly integrated into the readdcm_3d function.

{% raw %}

normalize[source]

normalize(t:Tensor)

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
normalize(im).mean(), normalize(im).std()

(tensor(4.8384e-08), tensor(1.))
{% endraw %} {% raw %} 
im = readdcm_3d('/media/ScaleOut/prostata/data/dcm/A0042197734/T2/DICOM', return_scaled=True, return_normalized=True)
im.shape, show_image_3d(im, axis = 0, nrow = 10)

((31, 736, 736), None)
{% endraw %}

Rendering 3D Objects

Somtimes the mask is better viewed as 3D object. Rendering is implemented as described in this example: https://scikit-image.org/docs/dev/auto_examples/edges/plot_marching_cubes.html A faster, more felxible way might be using ipyvolume or vtk.

To implement 3D rendering for the mask, the TensorMask3D class needs to be expanded

{% raw %}

class TensorMask3D[source]

TensorMask3D(*args, **kwargs) :: TensorDicom3D

Base class for 3d Segmentation, inherits from TensorDicom3D

{% endraw %} {% raw %}
{% endraw %} {% raw %} 
im  = TensorMask3D.create('/media/ScaleOut/prostata/data/dcm/A0042197734/T2/Annotation/Annotation.nii.gz')       
im.render_3d(alpha = (0.15, 0.15))

/home/bressekk/anaconda3/envs/fastai-v2/lib/python3.7/site-packages/ipykernel_launcher.py:127: FutureWarning: marching_cubes_lewiner is deprecated in favor of marching_cubes. marching_cubes_lewiner will be removed in version 0.19
{% endraw %} {% raw %} 
im.calc_volume()

{'background': 3694497.25, 'total_mask_volume': 25414.79155157879, 'class 1': 12753.419921875, 'class 2': 12661.4892578125}
{% endraw %}